home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / basic / marqfit.zip / MARQFIT.BAS next >
BASIC Source File  |  1985-12-11  |  34KB  |  669 lines
  1. 1 ' PROGRAM MARQFIT.BAS -- a non-linear least squares fitting
  2. 2 ' routine using the marquardt algorithm
  3. 5 ' (c) 1985 W. Schreiner, M. Kramer, S. Krischer, & Y. Langsam
  4. 20  '-------------- initialize program settings ---------------
  5. 30  REV = 2.1 : PI = 3.14159265#
  6. 40  ON ERROR GOTO 19000 : RETCOD = 0 :KEY OFF
  7. 60  DEF SEG=64:POKE 23,(PEEK(23) OR 64): DEF SEG 'set caps lock
  8. 65  ' The values MXVAR%, MXOBS% and MDI% can be changed-
  9. 70  ' be sure to change the DIM sizes correspondingly
  10. 80  MXVAR% = 30:MXOBS% = 500 : MDI% = MXVAR% * (MXVAR% + 1) / 2
  11. 90  DIM X(500),DTA(500),WGT(500),FXP(500),FYP(500),XP(500)
  12. 100 DIM YP(500),C(465),A(465),B(30),DPARM(30),PARM(30),                                 KFIX(30),DERIV(30),G(30),GRAD(30)
  13. 110 DEF FNLGT(X) = LOG(X)/LOG(10)
  14. 120 ON KEY(1) GOSUB 10200 : ON KEY(2) GOSUB 10300
  15. 130 ON KEY(3) GOSUB 10310 : ON KEY(4) GOSUB 10400
  16. 140 KEY(1) ON : KEY(2) ON : KEY(3) ON : KEY(4) ON
  17. 160   FG=2:BG=0:FG1=14:BG1=0:FG2=4:BG2=0 ' color attributes
  18. 165 FOR I = 1 TO MXOBS% : WGT(I) = 1 : NEXT
  19. 190 ' ---------------plot routine initialization
  20. 200 PLOTSET% = 0 ' plot definition flag
  21. 210 XTOTAL.PIX = 639: YTOTAL.PIX = 199
  22. 220 ' (for lo_res use XTOTAL.PIX = 319)
  23. 230 ' XMIN.PIX is lower left hand xmin in pixel coords;
  24. 240 ' YMIN.PIX is lower left hand ymin in pixel coords;
  25. 250 ' XMIN and YMIN are the same but in USER coordinates
  26. 260 XMIN.PIX=CINT(.11*XTOTAL.PIX): XMAX.PIX=.94*XTOTAL.PIX
  27. 270 YMAX.PIX=CINT(.07*YTOTAL.PIX)-2: YMIN.PIX=YTOTAL.PIX-12
  28. 280 DELTAX.PIX = XMAX.PIX-XMIN.PIX+1
  29. 290 DELTAY.PIX = YMIN.PIX-YMAX.PIX+1
  30. 310 ' functions to convert X & Y in user coords --> pixel coords
  31. 320 DEF FNSX(N)=XMIN.PIX + CINT((N-XMIN)*DELTAX.PIX/(XMAX-XMIN))
  32. 330 DEF FNSY(N)=YMIN.PIX - CINT((N-YMIN)*DELTAY.PIX/(YMAX-YMIN))
  33. 335 ' ------------------------start up--------------------------
  34. 340 IPFLG = 0    ' pause-in-fitting flag
  35. 350 ISVFL =  - 1 ' data saved flag
  36. 360 CLS : LOCATE 10 : COLOR FG2,BG2 :                                               PRINT TAB(15) "MARQUARDT LEAST SQUARES - REV ";REV:
  37. 365 COLOR FG,BG : FOR I = 1 TO 1700: NEXT:LOCATE 20
  38. 370 INPUT "USE DATA ON DRIVE (ENTER LETTER):";DRV$:                                 DRV$ =LEFT$(DRV$,1) : IF DRV$ <> "" THEN DRV$=DRV$+":"
  39. 380  GOTO 420
  40. 390 ' -----------------------   menu  --------------------------
  41. 400 PRINT "PRESS ANY KEY TO CONTINUE":WHILE INKEY$ = "":WEND
  42. 410 RETCOD = 0 ' reset error flag
  43. 420 F1%=0:SCREEN 0 :CLS: GOSUB 10000: COLOR FG2,BG: LOCATE 1,20
  44. 425 PRINT "MENU OPTIONS:": COLOR FG,BG : PRINT:COLOR FG1,BG1
  45. 430 PRINT TAB(6) "GENERAL:";:COLOR FG,BG
  46. 435 PRINT TAB(21) "1- ENTER TITLE " 'line  800
  47. 440 PRINT  TAB(21) "2- PRINTER ";   'line  850
  48. 450 IF PFLG% = 0 THEN PRINT "ON/";:COLOR FG1,BG1:PRINT  "OFF":                      COLOR FG,BG ELSE COLOR FG1,BG1: PRINT "ON";: COLOR FG,BG :                      PRINT  "/OFF"
  49. 460  PRINT  TAB(21) "3- SOLVE "      'line  900
  50. 470  PRINT  TAB(21) "4- PLOT "       'line  15000
  51. 480  PRINT  TAB(21) "5- QUIT "       'line  1300
  52. 490  PRINT  : COLOR FG1,BG1:PRINT  "     ENTER DATA:";:                              COLOR FG,BG:PRINT TAB(21) "6- MANUAL" 'line 1400
  53. 500  PRINT  TAB(21) "7- UREAD"       'line 1800
  54. 510  PRINT
  55. 520  COLOR FG1,BG1:PRINT  "     MODIFY DATA:";:COLOR FG,BG:                          PRINT TAB(21) "8- EDIT"         'line 2000
  56. 530  PRINT  TAB(21) "9- SCALE"       'line  2200
  57. 540  PRINT  TAB(20) "10- ZERO"       'line  2400
  58. 550  PRINT:COLOR FG1,BG1:PRINT "     PARAMETERS:";:COLOR FG,BG:                      PRINT TAB(20) "11- ENTER/REVIEW/CHANGE"     'line 2500
  59. 570  PRINT  "                   12- FIX"         'line  2800
  60. 580  PRINT  "                   13- FREE"        'line  3000
  61. 590  PRINT
  62. 600  COLOR FG1,BG1:PRINT  "     DATA & PARAM: ";:COLOR FG,BG:                        PRINT "14- LIST "     ' line 3200
  63. 610  PRINT  "                   15- SAVE ON DISK" 'line 3400
  64. 620  PRINT  "                   16- READ FROM DISK" 'line 3600
  65. 630  COLOR 15,BG: INPUT; "ENTER CHOICE";KMND:COLOR FG,BG:PRINT
  66. 650  ON KMND+1 GOTO 420,800,850,900,15000,1300,1400,1800,2000,                       2200,2400,2500,2800,3000,3200,3400,3600
  67. 660  PRINT  "*** ERROR *** ";KMND;" INVALID COMMAND": GOTO 400
  68. 800 'enter a title for documentation ---------------------------
  69. 810 INPUT "ENTER TITLE: ";TITL$:LOCATE 25,1:PRINT TITL$;SPC(11)
  70. 820  GOTO 420
  71. 850 'toggle printer on/off -------------------------------------
  72. 860 PFLG% = 1 - PFLG%: GOTO 420
  73. 900 'solve the mrqdt non-linear lst sq fit problem -------------
  74. 920 ITER% = 0 : IF IPFLG <  > 0 THEN ITER% = NITER%
  75. 930 IPFLG = 0:ISVFL = 0 : CLS : GOSUB 10000
  76. 940 INPUT "HOW MANY ITERATIONS?  ";MXITER%
  77. 950 IF MXITER% < 0 THEN MXITER% = 0
  78. 960 GOSUB 6000 ' go fit !
  79. 970 IPFLG =  - 1 : NITER% = NITER% + ITER%
  80. 990 PRINT:PRINT: PRINT TITL$ :PRINT NITER%;" ITERATIONS": PRINT
  81. 1000 IF PFLG% = 1 THEN LPRINT:LPRINT TITL$                                            :LPRINT NITER%;" iterations "
  82. 1010 GOSUB 9000
  83. 1020 IF FLG%=1 THEN PRINT:PRINT "CHOLESKY NEGATIVE DIAGONAL--";                            "UNABLE TO  SOLVE WITH SUPPLIED INITIAL PARAMETERS"
  84. 1030 PRINT  : PRINT  "PRESS ANY KEY FOR FITTING STATISTICS"
  85. 1040 WHILE INKEY$ = "" : WEND
  86. 1050 GOSUB 1100 : PRINT : GOTO 400 'print statistics and return
  87. 1070 'print fit statistics -------------------------------------
  88. 1100 PRINT  : PRINT  "SOME FITTING STATISTICS:"
  89. 1110 PRINT "SIGMA= ";SIG;" R= ";R;: PRINT TAB( 41);                                                                    "WGT'D SUM OF SQ. = ";WSS
  90. 1120 CHI2 =  INT (10 * CHI2) / 10
  91. 1130 PRINT  "CHI SQUARED= ";CHI2;" / ";NF;" DEG OF FREEDOM"
  92. 1140 PRINT  "# OF CALLS TO SUM OF SQ=";NSSC;:                                             PRINT TAB( 41);"# OF DERIV CALLS=";NDC:                                         PRINT  "# INC IN LAMBDA=";INCR;" LAMBDA= ";LAMBDA
  93. 1150 IF PFLG% <> 1 THEN RETURN
  94. 1160 LPRINT : LPRINT  "SOME FITTING STATISTICS:"
  95. 1170 LPRINT  "SIGMA= ";SIG;" R= ";R;:                                                 LPRINT TAB( 41);"WGT'D SUM OF SQ. = ";WSS
  96. 1180 LPRINT  "CHI SQUARED= ";CHI2;" / ";NF;" DEG OF FREEDOM"
  97. 1190 LPRINT "# OF CALLS TO SUM OF SQ=";NSSC;:                                        LPRINT   TAB( 41);"# OF DERIV CALLS=";NDC:                                      LPRINT  "# INC IN LAMBDA=";INCR;" LAMBDA= ";LAMBDA
  98. 1200 RETURN
  99. 1300 'quit -----------------------------------------------------
  100. 1310  IF ISVFL <  > 0 THEN 1360
  101. 1320  INPUT "PRESENT DATA NOT SAVED. SAVE IT?? (Y/N)";T1$
  102. 1330  IF  LEFT$ (T1$,1) = "Y" THEN 3400 'go save
  103. 1340  IF  LEFT$ (T1$,1) <> "N" THEN 420
  104. 1360  DEF SEG=64:POKE 23,(PEEK(23) AND 191):DEF SEG 'clear caps
  105. 1370  CLS: END
  106. 1400  'manual data entry ---------------------------------------
  107. 1410  PRINT  "MANUAL DATA ENTRY - ": PRINT  "   (, TO EXIT)"
  108. 1412 ' set flags to indicate data points that currently exist
  109. 1415 FOR K=1 TO NOBS%:XP(K)=1:NEXT :                                                  FOR K = NOBS%+1 TO MXOBS%:XP(K)=0:NEXT
  110. 1420 PRINT  "POINT    X,MEASUREMENT"
  111. 1425 PRINT  "===== ====,==========="
  112. 1430 FOR I = 1 TO MXOBS%
  113. 1440    PRINT  I;: INPUT " ";T1$,T2$
  114. 1450    IF T1$ = "" THEN NOBS% = I - 1: GOTO 1480
  115. 1460    X(I) =  VAL (T1$):DTA(I) =  VAL (T2$) :XP(I) = 1
  116. 1470 NEXT
  117. 1480 INPUT "WEIGHTS-ENTER 'NO','STAT' OR 'EXPLICIT'?";T1$
  118. 1490 IF T1$ = "NO" THEN FOR I=1 TO NOBS%:WGT(I)=1:NEXT:GOTO 1580
  119. 1500 IF T1$ = "STAT" THEN                                                              FOR I = 1 TO NOBS%: WGT(I)=SQR (DTA(I)): NEXT : GOTO 1580
  120. 1510 IF T1$ <  > "EXPLICIT" GOTO 1480
  121. 1520 PRINT  : PRINT  "POINT     X      Y      WEIGHT":                                        PRINT  "=====   =====  =====    ======"
  122. 1530 FOR I = 1 TO NOBS%
  123. 1540     PRINT  I;"     ";X(I);"    ";DTA(I);"     ";WGT(I);
  124. 1550     INPUT "New Weight=";WGT$:IF WGT$="" THEN 1570
  125. 1560     WGT(I) = VAL(WGT$)
  126. 1570 NEXT
  127. 1580 NOBS% = 0 : PRINT  "NOW ORDERING & REVIEWING DATA"
  128. 1600 FOR I = 1 TO MXOBS%
  129. 1610    IF XP(I) = 0 GOTO 1650  'zero it if pt is to be deleted
  130. 1620    NOBS% = NOBS% + 1
  131. 1630    X(NOBS%) = X(I):DTA(NOBS%) = DTA(I):WGT(NOBS%) = WGT(I)
  132. 1640    IF I <  > NOBS% THEN XP(I) = 0
  133. 1650 NEXT
  134.